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We propose a simple scheme to construct composition-dependent interatomic potentials for 
multicomponent systems that when superposed onto the potentials for the pure elements can 
reproduce not only the heat of mixing of the solid solution in the entire concentration range but 
also the energetics of a wider range of configurations including intermetallic phases. We show 
that an expansion in cluster interactions provides a way to systematically increase the accuracy 
of the model, and that it is straightforward to generalise this procedure to multicomponent 
systems. Concentration-dependent interatomic potentials can be built upon almost any type 
of potential for the pure elements including embedded atom method (EAM), modified EAM, 
bond-order, and Stillinger- Weber type potentials. In general, composition-dependent TV-body 
terms in the total energy lead to explicit N + 1-body forces, which potentially renders them 
computationally expensive. We present an algorithm that overcomes this problem and that 
can speed up the calculation of the forces for composition-dependent pair potentials in such 
a way as to make them computationally comparable in efficiency and scaling behaviour to 
standard EAM potentials. We also discuss the implementation in Monte-Carlo simulations. 
Finally, we exemplarily review the composition-dependent EAM model for the Fe-Cr system 
[PRL 95,075702, (2005)]. 
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potentials; cluster interactions 



1. Introduction 

Twenty-five years ago, the Finnis-Sinclair many body potential the Embedded 
Atom Model of Daw and Baskes 0] , the Glue model of Ercolessi and Parrinello 0] , 
and the effective medium theory due to Puska, Nieminen and Norskov (J @] marked 
the birthday of modern atomic scale computational materials science, enabling 
computer simulations at the multimillion atom scale to become a routine in modern 
materials science research. This family of many body potentials share in common 
the fact that the expression for the total energy has non linear contributions of pair 
functions, removing in this way the limitations of the pair potential formulation to 
describe realistic elastic properties. 

Alloys and compounds, where the thermodynamic information is of relevance, is 
one of the main fields in which these potentials have been applied. In the early days 
of many body potentials the main alloy property fitted was the heat of solution of a 
single impurity [a], i.e. the dilute limit of the heat of formation (HOF) of the alloy. 
However, when these potentials are applied to concentrated alloys the predictions 
are usually uncontrolled; they work well for systems with a mixing enthalpy that is 
nearly symmetric and positive over the entire concentration range, as for example 
in the cases of Fe-Cu @, Q, or Au-Ni @, 0, Ej. 

Alloys which show a strong asymmetry or even a sign inversion in the HOF such 
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as Fe-Cr or Pd-Ni are beyond the scope of standard many body potential models, 
and there is not yet a unique methodology suitable for their description. Similar 
limitations apply with respect to systems with a negative HOF which feature in- 
termetallic phases. Frequently, such systems require different parametrisations for 
different phases, as in the case of Ni-Al with the B2 phase on one hand [HI], and 
the 7 and 7' phases on the other 

Two schemes have been developed to deal with these shortcomings in the case 
of Fe-Cr which displays an inversion in the HOF as a function of concentration, 
namely the composition-dependent embedded atom method (CD-EAM) |13j ] and 
the two-band model (2BM) [13]. For neither one of these schemes, it is obvious 
how it can be extended to more than two components. 

The objective of this paper is to develop a framework for constructing interatomic 
potential models for multicomponent alloys based on an expansion in clusters of 
increasing size that can be practically implemented and systematically improved. 
Our methodology allows to describe systems with arbitrary heat of mixing curves 
and includes intermetallic phases in a systematic and physically meaningful fashion. 
Thereby, we overcome the most important disadvantages of current alloy potential 
schemes and provide a framework for systems of arbitrary complexity. 

In our methodology the interatomic interactions are modified by composition- 
dependent functions. This introduces a dependence on the environment which is 
somewhat reminiscent of the bond-order potential (BOP) scheme developed by 
Abell and Tersoff [H-[l3]. In this formalism the attractive pair potential is scaled 



by a (usually) angular dependent function (the "bond-order") which describes the 
local structure. Thereby, it is possible to distinguish different lattice structures 
(face-centred cubic, body-centred cubic, cubic diamond etc.) and also to stabilise 
structures with low packing density such as diamond or zincblende lattices. (In fact, 
the BOP formalism has been successfully applied to model alloys such as Fe-Pt that 
feature intermetallic phases with different lattice structures [18l|). The composition- 
dependent interatomic potential (CDIP) scheme introduced in the present work 
and the BOP formalism thus both include explicit environment-dependent terms. 
However, in the CDIP approach this environment-dependence is used to distinguish 
different chemical motifs while in the BOP scheme it is used to identify different 
structural motifs. 

This paper is organised as follows: In Sect. 12.11 we introduce the basic termi- 
nology and present a systematic approach to fitting potentials for binary systems. 
Section [2.21 describes how by including higher order terms it is possible to fit e.g., 
intermetallic phases. In Sect. [3TT1 a, series expansion is developed which generalises 
the concepts introduced in the previous sections and which is used in Sect. 13.21 
to obtain explicit expressions for a ternary system. The efficient computation of 
forces is discussed in Sect. [5] and an optimal implementation in Monte-Carlo simu- 
lations is the subject of Sect.[5l Finally, as an example, the composition-dependent 
embedded atom method potential for Fe-Cr is reviewed in Sect. [H 



2. Binary Systems 
2.1. Pair Potentials 

For the sake of clarity of the following exposition, we assume EAM models through- 
out this paper. It is important to stress that the formalism to be developed hereafter 
can be applied to any potential model for the pure elements including modified em- 
bedded atom method (MEAM) [H,[2Q|], bond-order (IE-E3], and Stillinger- Weber 
type [2l|] potentials. 
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Consider a single-component system of atoms A, whose interactions are described 
by the EAM model, 

E A = Y,U A {-p l ) + \Y.Y. ( t )A ^ with Pi = J2p( r n)- (!) 

i i j j^i 

The first term in Eq. ([T]) contains the embedding function [/^(pj), which is a 
nonlinear function of the local electron density 'p i around atom i. It accounts for 
cohesion due to band formation in the solid state and is constructed to reproduce 
the equation of state of system A. The second term represents the remainder of the 
interaction energy. It can be interpreted as the effective screened Coulomb inter- 
action between pairs of ions in A. The EAM formalism can capture the energetics 
associated with density fluctuations in the lattice and has been successfully applied 
for modelling the formation of crystal defects such as vacancies, interstitials and 
their clusters. 

Consider now a binary system, where the pure phases are described by EAM 
potentials. It can be shown that the total energy expression for this type of poten- 
tials is invariant under certain scaling operations [22|]. This "effective pair format" 
can be used to rescale the two EAM potentials, e.g. such that at the equilibrium 
volume for a certain lattice the electron density is 1, to ensure their compatibility. 
One part of the total energy of the two-component system can be written as the 
superposition of the respective embedding terms and effective pair interactions: 

eo=y,u a [pf + pa(b) pf) + \ E E ^ fa) ( 2 ) 

ieA ieAjeA 

+ J2 Ub + Pb{a) Pt) + g E E $ B fa) ' 

i£B i€B j&B 

where 

Pf= E Arv)- (3) 
jesj^i 

Note that above we have not yet added any explicit A—B interactions. Equation ([2]) 
is a strict superposition of the interatomic potentials for the pure elements with 
the only caveat that the electron density of the A (B) species in the embedding 
function of a B (A) particle is scaled with a parameter /j-b(A) i n order to account for 
the different local electron densities. Thereby, two EAM models can be calibrated 
with respect to each other. More elaborate schemes are possible, e.g. one can treat 
HA and /is as free parameters. Here for the sake of simplicity, we restrict ourselves 
to normalised electron densities. 

Starting from a parametrisation for Eq, we now devise a practical scheme for 
systematically improving the interaction model. Let us denote the true many-body 
energy functional of the binary system by Ef. Our goal is to construct an inter- 
atomic potential model for the difference energy functional AE^ = E t - E Q . We 
begin with the two dilute limits. Consider a lattice of A particles and substitute 
the atom residing in the i-th site with a B atom. Let us now assume that A£(°) 
for this configuration can be satisfactorily represented by a pair potential between 
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the A — B pairs. In this limit AE^ can thus be written as 

A^°)(^-rich) = ^^ B (r ii ). (4) 

(There is only one sum in this expression since we are dealing with a configuration 
that contains only one B atom). A similar expression is obtained for the I?-rich 
limit 

AE(°XB-nch) = J2vMn j )- (5) 

Since we do not require the pair potential models for the two dilute limits to 
coincide with each other, an interpolation is needed which preserves the energetics 
of the impurities. The main objective of the present paper is to devise such an 
interpolation scheme. The simplest ansatz for such an expression is 

A ^ (0) = E E 4 v Um) + E E 4 v AB(rv) (6) 

i£Aj£B iGAjeB 

Above, xfj denotes the concentration of species S in the neighbourhood of an 
A — B pair residing on the i and j sites. Ideally, we require this quantity to be 
easy to calculate and to be insensitive to the local density and topology, in other 
words it should separate chemistry from structure. In any case, xfj has to represent 
an average over the neighbourhood of both centres i and j. Before we derive the 
expression for xfj, it is instructive to discuss the corresponding one-centre quantity 
xf. It describes the local concentration of species S around atom i. A simple way 
to determine xf is to choose a local density function o"(r,j) and then to evaluate 
the following expression 

x s = Eu&s,j^i) a ( r ij) = of 

which is indeed rather insensitive to the local geometry. This is most obvious in 
the dilute limits. The local concentration xf at the site of an impurity atom % is 
either (if S is the minority species) or 1 (if S is the majority species) independent 
of the local structure. This is, however, strictly true only for the impurity atom. 
For the other atoms in the system x^ varies between and 1 depending on the 
distance to the impurity atom. Also for these particles, atomic displacements may 
change the value of Xj. A total decoupling of chemistry and structure is therefore 
not possible. The optimal choice for cr(ry ) is the function that minimises the effect 
of local geometry on xf . Although it is possible to choose different a- functions for 
the different species, we do not expect the quality of the final potential to depend 
crucially on the choice of a(rij). In fact, we expect the best choice for cr(rjj) to be 
the simplest one. 

It is now straightforward to extend Eq. ([7|) to define the concentration xfj in 
the neighbourhood of a pair of atoms residing on sites i and j. To this end, we 
first define a quantity x s , to represent the concentration of the species S in the 



(7) 
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cutoff for atom 2 

l-0«5o+a5o) = n =0.786 



-(x° +x° ) = A =0 214 

cutoff for atom 1 

Figure 1. Schematic illustration of the connection between xf and two-centre concentrations xf^ and 
their computation in a binary alloy according to Eqs. JjJ and (JSJ . Here, the cutoff function <r(r) which 
appears in Eq. J7) is assumed to be a step function which is 1 for r < r c and zero otherwise. 
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x l - „ 



neighbourhood of atom i excluding atom j: 

X S _ ^{k€S,k^i,k^j) a ( r ik) _ of - 5(S,tj)(r{rijj 



l-a( rij )/af 
x s I l-<r(rij)/Wi 



1 - cr(rij)/ai 



tj=S 



u = s 



where tj denotes the type of atom i, and 5(ti,tj) is 1 if i j = tj and zero otherwise. 
Using this quantity, the two-centre concentration xfj can be defined as follows 

x fi = o ( x f(7) + x m) ( 9 ) 



Hence, the two-centre concentration of the species S about the atom pair 
is the average concentration in the two separate neighbourhoods of sites i and j 
excluding both of these atoms. This definition, which is illustrated in Fig. Q3 has 
the important advantage that the interpolation scheme introduced in Eq. ([6]) does 
not modify the interactions in the dilute limits, since xfj is strictly or 1 in the 
two limits irrespective of the local structure. Furthermore, it is straightforward to 
generalise Eq. Q to multi-centre concentrations. For example, in the next section, 
we will explicitly discuss the construction of interatomic potentials using three- 
centre concentrations. 

Let us now revisit Eq. ([6]). As mentioned earlier this is the simplest ansatz for 
AEW that can reproduce the energetics of both dilute limits. A more general 
expression is 

A ^ (0) = E E h *Bi4) y Unj) + E E h A B (4) v AB(nj), (10) 

i€AjeB ieAjeB 

where h^(x) and h% are nonlinear functions with the property h^(0) = hg(0) = 
and frf (1) = = 1. By fitting these functions to the energetics of the 

concentrated alloys, the quality of the interatomic potential model for the binary 
can be improved drastically. 

In principle, one can stop here and have an interatomic potential model, Eq + 
A£(°), that can reproduce the energetics of the dilute limits as well as the solid 
solution of the binary. It is, however, also possible to further refine the above model. 
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For this purpose, let us again define a difference energy functional 

= E t - E — Ae(°\ (11) 

and construct an interatomic potential model for the energy functional AE^\ 
Consider a lattice of A particles and substitute two atoms, say i and j, with B 
particles. Assume that AE^ for this configuration can be well represented by 
a potential model describing the interaction of the B-B pair with a lattice of A 
particles. In this limit we can express A£« as 

A£«(,4-rich) = V£ B (r ij ) + ^V B A BA (n,k), (12) 

k 

where is shorthand for the three sets of positions of the i, j and k atoms, i.e. 
{rj, r.,-, rfc}. In the same way we obtain for the S-rich limit 

AE^(B-nch) = vUn,) + E V AAB(n jk ). (13) 

k 

Note that AE^ has both a two-body and a three-body component and thus can 
be decomposed as follows 

^=AE^ l + AE^ let . (14) 

In the next section we discuss how to incorporate the three-body contribution into 
the interatomic potential model. For now, we only consider AE^ ir . Following the 
same line of arguments that lead to Eq. (|10p . we obtain the expression 

A 41l = E E h BB(4) V BB(n 3 ) + EE h AA(4) V AA(^), (15) 
ieB j£B i£Aj£A 

which reproduces the contributions of the pair terms in the two limits given by 
Eqs. ()12p and ()13|) . The two non-linear functions have to fulfil the conditions 

>»aa(0)=^b(0) = (16) 
h% A (l) = hi B {\) = 1. (17) 

By fitting the functions h^ A and h BB in the intermediate concentration range to 
the energetics of the concentrated alloy, one can obtain a further improvement for 
the interaction model for the binary system. 

2.2. Beyond Pair Potentials 

In this section, we show that the formalism introduced in the previous section can 
be extended to multi-body interaction potentials, which enables us to capture the 
energetics of a wider range of phases including ordered compounds. In the previous 
section, we outlined a scheme to construct composition-dependent pair potentials 
for the potential energy landscape E + AE^ + AE^ . It was also observed that a 
proper formulation of AE^' requires incorporation of explicit three-body terms. In 
this section we describe how to construct such composition-dependent multi-body 
potentials. 
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Figure 2. Schematic illustration of the computation of three-centre concentrations in a binary alloy using 
Eqs. (0 and 12H . Here, the cutoff function o(r) which appears in Eq. J7J is assumed to be a step function 
which is 1 for r < r c and zero otherwise. 

First, we require an interpolation scheme to connect the two limits of the three- 
body term AJS^p let in Eq. (|14p . The simplest ansatz for such an expression is 



AE, 



(i) 

triplet 



E E E + E E E z&KSuifofc). ( 18 ) 

i£B jeB keA i&A jeAkeB 



where xfj k denotes the concentration of species S in the neighbourhood of the 
triplet residing on sites i, j and k. In analogy with the derivation of the two-centre 
concentration Eq. Q, we start from the one-centre concentration xf and define 
the intermediate quantity xf(jk) that represents the concentration centred around 
atom i excluding atoms j and k 



r s 

x i(jk) 



E 



cr(ru) 



of - 6(S,tj)a(rij) - 5(S,t k )a(r ik ) 



„5 



^2 {¥i>¥jt¥k) <?{ni) <Ji - a{ rij ) - a(r ik ) 

1 - [5{S, tjWnj) + (5(5, t k )a{r lk )] /of 



l-[a(rij) + a(r ik )]/(Ti 

and now following the same line of arguments leading to Eq. 
three-centre concentration xf- k as follows 



•^ijk 



J i(jk) 



S , S 

x j(ik) + x k(ij) 



(19) 
(20) 

we define the 
(21) 



A graphical illustration of the computation of this quantity is given in Fig. [2j The 
three-centre concentration of the species S about the triplet (i,j,k) is the average 
concentration (excluding the triplet) in three separate neighbourhoods, each of 
which is centred at one of the atoms in the triplet. Thanks to this definition xfj k is 
strictly or 1 in the two dilute limits described in Eqs. (112j) and ([13]), irrespective 
of the local structure. Hence, the interpolation scheme in Eq. (I18j) does not alter 
the interactions in Eqs. (112j) and (j 13H . Again, as in Eq. (jlOp . we can improve the 
simple interpolation scheme in Eq. (|18|) 

A-^triplet = E E E ^BBA( x tjkWBBA( r ijk) + E E E ^AAB{ x ijk)^AAB{ r ijk) > 
ieBjeBkGA ieA jeAkeB 

(22) 

where h BBA and h AAB are non-linear functions that can be fitted to the energetics 
of the concentrated alloys with the boundary conditions 



hi BA {Q) = ^ab(O) = and h BBA (l) 



Kab(1) 



1. 



(23) 
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• -embedded 
cluster of order 2 
(•O) 



• -embedded 
cluster of order 3 
(•0») 



• -embedded 
cluster of order 4 
(•O 2 *) 

(C-^ ,C2,C 3 ) 



Figure 3. Schematic illustration of S-embedded coloured clusters of orders 2, 3, and 4 in a ternary alloy. 
The shaded region indicates the cutoff range around the central atom marked by an asterisk. 



Following this scheme composition-dependent cluster interactions of arbitrary order 
can be included in the interatomic potential model. To summarise, to incorporate 
cluster interactions of order n, two cluster potentials are constructed, one for the 
configuration where the cluster is embedded in the A lattice and one for the config- 
uration where the cluster is embedded in the B lattice. Subsequently these limits 
are interpolated using the n-centre concentrations. In the next section, we review 
this strategy in detail to show that a systematic series expansion in composition- 
dependent cluster interactions is possible for general multicomponent systems. 



3. Multicomponent Systems 

3.1. Series Expansion in Embedded Cluster Interactions 

In the first sections of this paper, we have shown how to practically construct 
interatomic potentials for binary systems. First, mixed interatomic pair and triplet 
potentials are generated for the dilute limits which are subsequently extended to 
arbitrary concentrations by fitting interpolation functions that depend on the local 
concentration about the atomic pairs and triplets. The choice of specific potentials 
and dilute configurations was mainly driven by physical intuition. In this section we 
show that this procedure can be formalised and generalised to arbitrarily complex 
systems with more than two components. 

Consider an n-component mixture of N particles that are distinguishable only 
through their species. Assign a unique colour to each of the species: {Ci, . . . ,C n }. 
We define a colour cluster of order m to be a set of m particles with a specific 
colour combination. We use the occupation number formalism to identify colour 

schemes, i.e. (c^ 1 , . . . , C^"^ , where ki is the number of particles in the cluster with 

colour Ci, and Y2i ^« = m - For example, a cluster of order 3 consisting of one parti- 
cle with the colour C\ and two particles with the colour C3, is denoted by (Ci,C|). 
Furthermore, we define an <S-embedded colour cluster of order m to be a set of 
m coloured particles embedded in a pure matrix of species S. Three examples of 
such 5-embedded coloured clusters are shown in Fig. [3l The key idea is that the 
potential energy landscape of an alloy can be expanded in the basis set of elemen- 
tary interaction potentials each of which is constructed to reproduce the energetics 
of a particular embedded colour cluster. The order of an interaction element in 
the series is determined by the order of the corresponding colour cluster. By pro- 
gressively including higher order colour cluster interactions, one can systematically 
increase the accuracy of the model. 

To recapitulate, we expand the potential energy landscape of multicompo- 
nent systems in the basis set of colour cluster interatomic potential functions 
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j(x* m +x* a) ) = l 4 =0.500 



2 ( x im +x 2tu ) ~~ 14 —0.214 



!(*•,+*•» )=n =0 - 286 

Figure 4. Schematic illustration of the connection between xf and two-centre concentrations xf^ and 
their computation in a ternary alloy according to Eqs. J7J and (5). Here, the cutoff function <r(r) which 
appears in Eq. (0 is assumed to be a step function which is 1 for r < r c and zero otherwise. 

c «,„ ({r}), where {r} is the real-space configuration of the respective clus- 
ter. The expansion coefficient for each basis function is the interpolation function 
h^ kl cfc „(a;' 5 ), where x s is the local concentration of the species S in the neigh- 
bourhood of the cluster. One of the innovations in this work is a simple and com- 
putationally expeditious way to determine x s which is illustrated for the case of 
a ternary alloy in Fig. [H Formally the total energy expression for an alloy of n 
components and N particles can be written as 

£ = £ o + EE-EE h &..s* .rs- (M), (24) 

m ki k„ S 

V v ' 

where the first sum is over the order of the cluster potentials and the subsequent 
sums are over all distinguishable colour combinations of m-size clusters. Each term 
in the above expansion can be evaluated as follows 

N N 

^ l e kl ...Cn n ( x<S )^cf 1 ...Cn n ({ r }) = ■ ■ ■ ^ ^c^.-.c*™ ( x ?i---« m )^cf i ...c* n ( r *i---* m )- 

«1=1 im = l 
V v ' 

{i 1 ...i m }e{c1 1 ...c^} 

(25) 

The sums in Eq. (|25j) are over all possible m-size atom clusters {i\ . . . i m } in the 

system with the colour scheme 1 , . . . , C^ n J . 

The main advantage of this scheme is that the basis functions can be constructed 
sequentially and independent of the interpolation functions. The lower order terms 
can be constructed with no knowledge of the higher order terms and therefore need 
not be reparametrised when higher order cluster potentials are constructed. The 
higher order terms in the expansion become progressively smaller. Furthermore, 
addition of new terms in the series expansion is not likely to introduce unphysical 
behaviour, a problem that plagues most fitting schemes for interatomic potentials. 

3.2. Explicit expressions for ternary alloys 

In this section we illustrate the formal discussion in the previous section by con- 
structing an expansion in embedded pair and triplet potentials for a ternary system. 
For simplicity we assume the pure elements are described by EAM models. The 
extension to larger number of components and higher order cluster potentials will 
be obvious. We consider a system of three components A, B and C, and assume 
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that three composition-dependent pair potentials for the binary systems A — B, 
A — C and B — C have already been constructed. Explicitly, the A — B interaction 
is given by the following expression 

EZb = E u * tpt + Pf) + \ E E (^04)^) + fcii(xg)V^(r -)) 

ieA ieAjeA 

+ E ^ + PB(A) P?) + j E E 01b(^)<^M + ^(4)V^ B (r^)) 

ieB ieB jeB 

+ E E + fcfB(xg)V&(r«)) • (26) 

ieA jeB 

By now the notation above should be familiar. The interaction potentials for the 
two other pairs can be written in analogous fashion. 

Now, we can spell out the expansion in embedded pair potentials for the ternary 
A-B-C 

ET B -c = E ^ (Pt + PMB) pf + Pa { c) Pf) (27) 

ieA 

+ E Ub + Pb(A) pf + PB(C) Pi) 

ieB 

+ E U ° (P? + PC( A ) ~Pi + PC(B) pf) 

iec 

+ \ E E MaO**})^) + *ii(*g)Vft(ry) + /&(zg)V&M 

+ \ E E [^b(^)<M%) + ^B(^)^B B (rv) + ^b(4)^bb(^)] 
ieB jeB 

+ \ E E [^c(*g)te(ry) + ^c(^)^co(r«) + ^c(^)^oc(r«)] 
iecjec 

+ E E MB(^)V^(r -) + /»2b(^)V&(^-) + h% B (x%)V$ B ( rij )] 
ieA jeB 

+ E E Wc(4)4N) + h B AC (x^ c (r l3 ) + ^(sgjV&fa)] 
ieA jec 

+ E E [^(4)^(^0 + ^c(^)V-#c(^-) + hZcixfyvgdrv)] . 
ieBjeC 

The only unknowns in the above equation are V^ B (rij), V AC (rij), V B( -,(rij), 
h AB( x ?j)> h Ac( x fj) and h Bc( x tj)- The Potentials V£,(r -), V&faj) and V^faj) 
describe the interaction between pairs of unlike species embedded in pure lattices 
of the third species of the ternary. In analogy with the previous section, it is reason- 
able to expect that we can construct these potentials separately in their respective 
dilute limits and subsequently fit the interpolation functions h < ^ B (xfj), h AC (x B ), 
h BC (xfj) to the energetics of the concentrated ternary alloys. However, when the 
number of species increases certain complications can arise that are not present in 
the binaries. This is well illustrated in the situation above. We now show that it 
is in fact not possible to separately construct the three pair potentials V AB (rij), 
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Figure 5. Schematic illustration of two and three-centre concentrations for a ternary alloy in the dilute 
limit. Note that the two-centre concentrations Xij in the dilute limit in a binary alloy are either one or 
zero. In contrast, in the case of a ternary alloy the two-centre concentrations in the same limit can be 
non-zero. The three-centre concentrations, however, are again either one or zero. 



V^ G (rij) and Vg C (rij) described above. 

To this end, consider a pure lattice of N particles of e.g., C species. Substitute 
two nearest neighbour particles in this lattice with an A particle and a B particle 
respectively. The ternary energy Eq. (|27p for a C-rich configuration containing one 
A — B pair on the sites i and j respectively becomes 

E p ^ B _ c (C-nch) = E 

+ \ E E ( h CC^c(r kl ) + h£ c (xi)v£ c (r kl ) + h^ c (x B kl )Vi c {r kl )) 
keciec 

+ VA B {n 3 ) + {h B AC {xf k )V% c (r lk ) + h c AC {xg)v2 c {r ik )) . 

k&C 

+ E { h Bc{4k) V Bc(r jk ) + h B ! c(xf k )Vg c (r jk )) , 
fceC 

where for the sake of clarity we have replaced the three embedding terms in Eq. ()27|) 
by Eq. Observe that all three unknown potentials V AB (rij), V B c (ri k ) and Vsc( r jk) 
as well as their corresponding interpolation functions appear in Eq. (|28p . This is in 
contrast to the binary case, e.g. Eqs. (|3J), ©, (fT2|) and (fT3|) . where the potentials 
for the two dilute limits can be constructed independently of each other. This is 
because the two-centre concentrations in the dilute limit in a binary alloy are either 
one or zero. In contrast, in the case of a ternary alloy the two-centre concentrations 
in the same limit can be non-zero (see Fig. [5]). 

A straightforward solution to the above problem is to fit all three pair potentials 
simultaneously. A closer look at Eq. (|28p . however, suggests a simpler solution. Let 
us examine the interpolation functions h B c (xf k ) and h BC (x^ k ). Note that since we 
are dealing here with an A — B cluster in a C-rich system xf k and x^ k are close 
to zero. Remembering the boundary conditions on the interpolation functions, i.e. 
h(l) = 1 and h(0) = 0, we conclude that the contributions of the V B c (rij) and 
V BC (rij) potentials to the energetics of an A — B pair embedded in a C lattice are 
small. In fact, we can diminish the contribution of these potentials to Eq. (|28p by 
enforcing the interpolation functions to be for x < xth, where xth is the largest 
concentration of B or A particles found about any pair in the system. In this way, 
one can generally separate the construction of cluster potentials when they overlap 
in the dilute configurations. 

The problem of potential overlap in the dilute limit discussed above should not 
be neglected. On the other hand it is quite benign and — as shown above — can be 
handled easily. Furthermore, more often than not, even for complex clusters and 
many components, there is no overlap. We illustrate this point by considering the 
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simplest expansion in triplet cluster potentials for the ternary above: 

E^b-c = E E E h A ABo{4k) vAcinjk) (28) 
ieA jeB keC 

+ h-ABc( x Fjk) VABc( r ijk) + h^Bc( x ?jk) ^ABc( r ijk)- 

Now consider again the same C lattice as above, where an A — B pair has been 
embedded at the sites i and j. The triplet energy becomes 

£« c (C-rich) = Y, yABc(r ijk ). (29) 
kec 

Since we have only contributions from ^ABc^ijfc) ^ or ^ ne these configurations, 
we can construct these potentials separately from each other and independent of 
the interpolation functions. This is because in the dilute limit the three-centre 
concentrations are again either one or zero (see Fig. [5]) . 



4. Implementation of Forces in Molecular dynamics 

Next to accuracy, the most important quality of an interatomic potential model 
is its computational efficiency when implemented into atomistic simulation codes. 
Due to the unconventional form of the interatomic potentials described in this 
work, it is important to discuss the efficient implementation of forces for molecular- 
dynamics simulations. We will see below that the straightforward derivation of the 
forces for composition-dependent pair potentials leads to explicit 3-body forces. 
In fact in general, composition-dependent TV-body potentials lead to explicit N + 
1-body forces. Below we present an algorithm that considerably speeds up the 
calculation of forces for composition-dependent iV-body potentials, making them 
comparable in efficiency to the corresponding iV-body regular potentials. In the 
following, for the sake of clarity we limit our discussion to pair potentials. The 
extension to cluster potentials of higher order is straightforward. 

For reference, let us first consider a conventional mixed pair potential energy 
expression for a binary system, 

^pp = EE y M- ( 3 °) 

ieA jeB 

Within this model the force on a particle k of type A is calculated as follows 

f? = £n^ (3D 

OT k jeB r kj 

Let us now consider a typical composition-dependent pair potential model for the 
same binary system, 

^cdp P = EE /i (4)^)' ( 32 ) 

i£A j£B 

where xf, is the two-centre concentration of the species A about the pair. 
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Now the force on particle k of type A can be written 



^ = E vtawxft) + EE%w4)^ + 5fJ ■ (33) 



for which after some algebra we obtain 

All the quantities above have already been defined in Eqs. (JSj) and ([9]). The second 
term in Eq. (|33p contains contributions from two particles i and j to the forces on 
particle k. Hence composition-dependent pair potentials lead to explicit three-body 
forces, which usually implies significantly more expensive to calculations. However, 
we will now show that in the case of expressions such as Eq. (j33[) one can regroup 
the terms in such a way as to speed up the calculation of forces drastically. To this 
end, let us introduce a per-atom quantity that for an atom of type A reads 

Mg A = £ V( rij )h'(4) ^^^ , (35) 
and for an atom of type B 

Mg^Yv^i ^^^ . (36) 

Substituting Mf into Eq. ([33]) we obtain 

^BE = £ V\r kj )h(4 3 ) + \ E Mfo\r h ^. (37) 

Similar derivation for the force on a particle k of type B leads to the expression 



^-SP = E + 5 E *Mr«)£. (38) 



Each quantity in the above force expressions can be calculated separately via pair- 
wise summations. This allows for a very efficient three-step algorithm for the cal- 
culation of forces: (i) compute and store the local partial densities of for every 
atom, (ii) compute and store the quantities Mf for every atom, and (Hi) compute 
the forces according to the Eqs. ([37|) and (i38j) . This method leads to computational 
efficiency comparable to standard EAM models. 



5. Linearised Models for efficient Monte-Carlo simulations 



Molecular dynamics simulations are limited when it comes to modelling phenomena 
such as precipitation, surface and grain boundary segregation, or ordering in alloys. 
Monte-Carlo (MC) methods, however, are ideally suited for such applications. The 
most common techniques are based on so-called swap trial moves, in which the 
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chemical identity of a random particle is changed. The resulting change in potential 
energy, AE, is used to decide whether the swap is accepted or rejected. 

The main task in an MC simulation is therefore to calculate the change in po- 
tential energy induced by swapping the type of a single atom. For short-range 
potentials this can be done very efficiently, since the type exchange only affects the 
atoms in the neighbourhood of the type swap. In the framework of the standard 
EAM model the situation is as follows: Changing the species of one atom directly 
affects (1) its embedding energy, (2) its pair-wise interactions with neighbouring 
atoms, and (3) indirectly changes the electron density at neighbouring atoms and 
therefore their embedding energies. All these quantities need to be recalculated by 
visiting the atoms affected by the type swap. 

In the case of composition-dependent models the situation turns out to be more 
laborious. To illustrate this let us again consider a typical composition-dependent 
pair potential model for a binary system: 

£cdpp = E5>(4) y M> ( 39 ) 
ieA jeB 

where xf- is the two-centre concentration of the species A about the pair 

where the is the local concentration A about the atom i excluding atom j. 
From Eq. ([9]) we observed that to a good approximation Xj. Therefore, 

for the qualitative discussion below we replace x^ by X{. In the energy expres- 
sion Eq. (f39j) . the site energy Ei of an atom i does not only depend on the local 
concentration Xj, but also on the concentrations Xj of all its neighbours j. This 
has a dreadful impact on the efficiency of the energy calculation. Changing the 
chemical identity of some atom i alters the local concentrations Xj of all its direct 
neighbours j, which in turn affects the mixed interaction of all atoms j with all of 
their respective neighbour atoms k. All of these have to be re-evaluated to compute 
the total change in energy induced by the single swap operation. The interaction 
radius that has to be considered is therefore twice as large as the cutoff radius of 
the underlying EAM potential, which increases the computational costs by at least 
one order of magnitude. 

This issue can be resolved quite effectively if we linearise the interpolation func- 
tion h(Xjj) as follows 

Hxi j ) = ±(h(x? {j) ) + h(xf {{) )). (41) 

Within the new linearised formulation, although a single pair interaction between 
two atoms j and k still depends on the concentration at both sites, the site energy 
can be recast in a form that is independent of the concentrations on the neighbour- 
ing sites. As a result, the site energy of atom k is no longer affected by changing 
the type of an atom i that is farther away than one cutoff radius. Note that lin- 
earisation can be done for interpolation functions of any n-centre concentrations. 
All composition-dependent models independent of cluster size can therefore be lin- 
earised. We have discussed the linearised model and its implementation for MD 
and MC at length in a recent publication [23]. 
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6. A practical example 

To provide a practical illustration of the concepts developed in this paper, we now 
revisit the composition-dependent EAM potential for Fe-Cr [3], which has already 
been successfully applied in a number of cases [H, [HI]. 

6.1. Application of composition- dependent embedded atom method to Fe-Cr 

Iron alloys are materials with numerous technological applications. In particular 
Fe-Cr alloys are at the basis of ferritic stainless steels. It has been recently shown 
[1^] that the Fe-Cr alloy in the ferromagnetic phase has an anomaly in the heat of 
formation which shows a change in sign going from negative to positive at about 
10% Cr and leads to the coexistence of intermetallic phase [27( and segregation in 
the same alloy. This c omp lexity results from a "magnetic frustration" of the Cr 
atoms in the Fe matrix [28j] which leads to an effectively repulsive Cr-Cr interaction. 
Capturing this complexity with an empirical potential model has been an active 
subject of research in recent years. 
To model this system, Caro and coworkers used the following ansatz 



£Fe-Cr = £ U Fc (pf e + p?") + \ £ ^ M ( 42 ) 

ieFe igFejGFe 

ieCr ieCrjeO 

+ £E /1 ( x -^p) 

where we used the same notation as in the earlier sections. The partial electron 
densities pf follow the same definition as in Eq. ([3]). Furthermore, the local con- 
centration variable X{ in Eq. (I42p is defined as 

^ = Grl Fe ' ( 43 ) 
Pi + Pi 

The two densities p Fe (rij) and p {rij) are normalised such that at the equilibrium 
lattice constant of each pure lattice, the respective partial electron density is 1. In 
this way the two EAM models for the pure elements are made compatible with 
each other. 

Equation (I42p looks quite similar to the composition-dependent pair potential 
energy expressions discussed in Sect. 12.11 There are, however, three essential dif- 
ferences: (i) There is only one mixed pair potential V m i xe d(rij) as opposed to two 
in Sect. 12.11 (one for each limit), (ii) There is no boundary conditions on the in- 
terpolation function h(x) at x = and x = 1. (iii) The local concentration about 
the pair is just the average of the one-centre concentrations about the two 
sites, and not the two-centre concentration as defined in Eq. ([9j). Of course, at no 
extra cost the more rigorous definition in Eq. ([9]) is a better choice for the measure 
of local concentration about a pair of atoms. On the other hand, Eq. ([8]) shows 
that the one-centre concentration above is only a perturbation away from the more 
accurate quantity. 

The Fe-Cr CD-EAM model was the pioneering work that has inspired the current 
paper. Here, we have tried to give a more rigorous foundation to the CD-EAM 



January 31, 2012 2:3 Philosophical Magazine prim 



3386 



B. Sadigh, P. Erhart, A. Stukowski, and A. Caro 



3 6 ^ 

<D 

E 

'% 4 
o 

1 2 
_o 

O 



linearized CD-EAM (Fe 50% Cr) 
original CD-EAM (Fe 50% Cr) 
Standard EAM (pure Fe) 

-« * 



27 64 
Number of processors 



512 



Figure 6. Comparison of the computation times for the CD-EAM models and the standard EAM model 
in a parallel molecular dynamics simulation. The benchmark simulation consists of a body-centred cubic 
crystal at 300 K with 16,000 atoms per processor. 

model. In fact, we can strictly argue that CD-EAM is a simplified version of the 
current formalism. It works very well for the Fe-Cr system since the two elements 
are similar in size and chemical nature. It is therefore reasonable to make the 
approximation that functional forms of the mixed pair potentials describing the 
two dilute limits are the same. 

Let us illustrate the last statement with the example of Lennard-Jones (LJ) 
potentials. These potentials are determined by two parameters: a and e; the first 
parameter specifies the position of the minimum of the potential or in other words 
the particle size, and the second parameter specifies the interaction strength. A 
mixture of two types of LJ particles with no size mismatch (same a) but different 
cohesive energies can be described by the same potential that is merely scaled 
differently for the two particles. Extending this analogy to the case of the Fe-Cr 
system we can see why only one mixed potential can be enough. However, it is 
important to realise now that when only one potential is used, the functions h{x) 
provide the interaction strength, which in the case of Fe-Cr is positive in one 
dilute limit and negative in the other. Hence no boundary conditions exist at the 
two concentrations x = and x = 1. 

In the original CD-EAM model, there was a further simplification. The mixed 
potentials V m \ x {rij) was never fitted. In fact it was taken as the average of the effec- 
tive EAM pairwise interactions of the pure elements at their respective equilibrium 
volumes 



1 



= ^ (0Fe(ry) + 2U Fe (pl e )p Fe ( rij ) + Mrij) + 2U C r(p$ I )p Cl "(ry)) , (44) 



where Pq is the electron density at the equilibrium lattice constant for the species 
S. Only the function h(x) was fitted to the heat of mixing of the solid solution. 
The success of this model in spite of all the simplifications is a telltale of the power 
of this methodology. 



6.2. Molecular dynamics and Monte Carlo performance 

In Sect. H] we presented an algorithm for calculating forces within the composition- 
dependent interatomic potential models which brings their efficiency on par with 
the standard EAM scheme. This was first discussed in a recent publication by the 
present authors [23| , where this algorithm was implemented for the Fe-Cr CD-EAM 
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Figure 7. Comparison of the timing in a MC simulation of a Fe-Cr alloy at 50% composition. The 
simulation cell contained 1024 atoms. 

model in the popular massively-parallel MD code LAMMPS [l^j.To benchmark its 
performance, we carried out MD simulations of a body-centred cubic (BCC) crystal 
at 300 K using periodic boundary conditions. For the CD-EAM case we considered 
a random alloy with 50% Cr. For the standard EAM case, the sample contained 
only Fe. Simulations were run on 1, 8, 27, 64, and 512 processors with 16,000 
atoms per processor (weak scaling). The results for the CD-EAM routines and 
the LAMMPS standard EAM routine are displayed in Fig. [6l In this figure, the 
original CD-EAM model as well as its linearised version are displayed. We see that 
the two versions are between 60% (linearised model) to 70% (original model) slower 
than the standard EAM. This is a small price to pay considering the fact that the 
CD-EAM expression actually contains explicit three-body forces. 

In our recent publication [23] we also studied the Monte Carlo performance of 
composition-dependent interatomic potentials focusing on the comparison of the 
original and the linearised CD-EAM model. The performance gain due to the lin- 
earised formulation is illustrated in Fig. [7] which compares the timing of the lin- 
earised and original CD-EAM models in a serial MC simulation for a random Fe-Cr 
alloy at 50% composition. We find that the linearized CD-EAM model is twelve 
times faster than the original formulation. This is an impressive performance gain, 
which clearly advocates for linearised composition-dependent interatomic poten- 
tials. 



7. Conclusions 

The present work has come about in response to a need for a practical scheme 
for fitting interatomic potential models for multicomponent alloys. At this point 
of time, when faced with the task of modelling the chemistry of e.g. a ternary 
alloy, one is overwhelmed with the complexity of the problem. In this paper, we 
have presented a systematic methodology for the construction of alloy potentials, 
starting from pre-existing potentials for the constituent elements. The formalism 
represents a generalisation of the approach employed by one of the authors for 
the Fe-Cr system [3]. We have shown that this formalism naturally extends to 
treating multicomponent systems. The main idea of the approach is to describe 
the energetics of dilute concentrations of solute atoms in the pure host in terms 
of pair and higher-order cluster interactions (see Figs. [3] and E]). These interaction 
functions are then used as a basis set for expanding the potential energy of the alloy 
in the entire concentration range. To describe the energetics of the concentrated 
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Figure 8. Several examples for clusters used to construct higher-order interaction terms which can be 
extracted from the configuration shown on the left. 



alloys, the contributions of the basis functions are weighted by interpolation func- 
tions expressed in terms of local concentration variables. One of the innovations in 
this work is a novel measure of local composition around individual atoms in the 
system. This introduces an explicit dependence on the chemical environment. In 
this sense the composition-dependent interatomic potential scheme is reminiscent 
of the bond-order potential scheme developed by Abell and Tersoff HI- 17] which 
employs a measure of the bond-order to distinguish between different structural 
motifs. 

The main advantage of the framework presented here is that the basis functions 
can be constructed sequentially and independent of the interpolation functions, 
leading to a scheme that can be practically implemented and systematically im- 
proved upon. The lower order terms can be constructed with no knowledge of the 
higher order terms and therefore need not be reparametrised when higher order 
cluster potentials are constructed. The higher order terms in the expansion become 
progressively smaller. In this way the model can be made step by step, starting 
from the lowest order cluster potentials. Furthermore addition of new terms in the 
series expansion is not likely to introduce unphysical behaviour, a problem that 
plagues most fitting schemes for interatomic potentials. 

The practical determination of the basis functions and the interpolation functions 
proceeds by fitting to first-principles data. The expansion in cluster interactions 
may be reminiscent of the celebrated "cluster expansion" technique (3(J that has 
been used extensively during the past few decades to model the thermodynamics of 
multicomponent alloys from first principles. But it is important to note here that 
the methodology presented in this paper has no relation to the cluster expansion 
technique. The latter reduces the continuous phase space of e.g., a binary alloy 
onto the discrete configuration space of the corresponding Ising model. There is 
only one number associated with each cluster configuration, namely the the free 
energy of that cluster. The so-called "effective cluster interactions" (ECIs) are 
usually obtained via an optimisation process from all the cluster free energies. A 
procedure of the sort proposed in this paper is not possible, since there is not direct 
link between any single cluster free energy and an ECI. In contrast, when fitting 
e.g. a VABi^ij) interaction potential, a solute inclusion not only changes the total 
energy of the system, it causes forces in the system and modifies the force constants 
of the host, all of which can be used to construct a continuous pair potential. 

Composition-dependent interatomic potentials are constructed by incorporating 
pair, triplet and higher-order cluster interactions that describe the energetics of 
clusters embedded in a pure host with a specific underlying lattice. One may now 
wonder, with this approach, could a potential be expected to handle systems which 
change lattice type as a function of concentration? For instance the Ni-Al phase 
diagram contains phases with BCC-based crystal structures, while the pure metals 
are face-centred cubic (FCC). Following the approach described above, the basis 
functions are parametrised in terms of solute cluster energies in the constituent 
FCC structures. How can one then expect to provide a reasonable model for the 
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BCC-based NiAl phase? The answer lies in the interpolation functions. They are 
fitted to the energetics of the ordered and disordered compounds along the con- 
centration range with arbitrary crystal structures. 
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